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ABSTRACT 

We perform a series of simulations of evolving star clusters using AMUSE (the Astrophysical 
Multipurpose Software Environment), a new community-based multi-physics simulation package, and 
compare our results to existing work. These simulations model a star cluster beginning with a King 
model distribution and a selection of power-law initial mass functions, and contain a tidal cut-off. They 
are evolved using collisional stellar dynamics and include mass loss due to stellar evolution. After 
determining that the differences between AMUSE results and prior publications are understood, 
we explored the variation in cluster lifetimes due to the random realization noise introduced by 
transforming a King model to specific initial conditions. This random realization noise can affect 
the lifetime of a simulated star cluster by up to 30%. Two modes of star cluster dissolution were 
identified: a mass evolution curve that contains a run-away cluster dissolution with a sudden loss of 
mass, and a dissolution mode that does not contain this feature. We refer to these dissolution modes 
as "dynamical" and "relaxation" dominated respectively. For Salpeter-like initial mass functions, we 
determined the boundary between these two modes in terms of the dynamical and relaxation time 
scales. 



1. INTRODUCTION 

Star clusters are natural laboratories for many astro- 
physical processes. In the simplest description, cluster 
stars may be thought of as being (almost) coeval point 
masses — an N-body system — and their motion traces 
their mutual gravitation and the possible influence of an 
external galactic tidal field. In more complex situations, 
stars evolve, gas may accrete into the cluster, new stars 
may form out of that gas, and the gas may be expelled 
from the cluster quickly by supernovae or more slowly by 
radiation pressure and stellar winds. A typical cluster is 
subject to several long-term mass-loss processes, includ- 
ing losses due to stellar evolution and tidal stripping of 
the outermost stars. These processes compete with re- 
laxation processes to define the equilibrium state of the 
cluster. 

Setting aside the complexities of intracluster gas, sim- 
ple models combining a few basic physical processes — 
stellar dynamics, stellar evolution, and tidal effects — 
have proved very useful in the study of star clusters. 
These simulations combine differing treatments of multi- 
ple physical processes, and must be carefully calibrated 
to ensure their reliability. 



Chernoff & Weinberg 



(1990) 



(noted as "CW" throughout this paper) combined a sim- 
ple stellar evolution prescription with Fokker-Planck sim- 
ulations of stellar dynamics and a highly idealized tidal 
field to produce a seminal "ba seline" set o f cluster sim- 
ulations, starting from King (King 1966) initial mod- 



els. This survey, together with subsequent studies by 
Fukushige fc Heggiel dl995t , |Aarseth fe Reggie] ( |1998| ), 
and [T'akahashi ^FToiTegies Zwart| ( |2000p (abbreviated 
as '"i'FZ" for the remainder of this paper), using other 
formulations of stellar evolution and both N-body and 
Fokker-Planck treatments of stellar dynamics, have re- 
sulted in comparative catalogs of parameter space that 
now serve as tests of any new code. 

Part of the purpose of this paper is to validate parts of 
AMUSE [j] — the Astrophysical Multipurpose Simulation 
Environment — against known results and then to show 
new applications of the framework to stellar cluster dy- 
namics. AMUSE is a new software framework designed 
for simulations of dense stellar syst ems, inspired by the 
earl ier MUSE project described by |Harfst et"aL] ( 2008 1 
and Portegies Zwart et al. (2009). A detailed technical 
acco unt of AMUSE is beyond the scope of this articl e 
(see |McMillan et al.||2011| | Portegies Zwart et al.l|2012[) . 
A summary is presented in jJ5]to provide the reader with 
some context on the software used. 

We set out to test AMUSE against known results, 
but found that comparing different simulations at any 
meaningful level of precision is a non-trivial task. In this 
paper, we employ an N-body stellar dynamics code, sev- 
eral stellar evolution codes, and a simple tidal stripping 
algorithm as the three basic simulation components, and 
compare AMUSE with the results of TPZ. 

This line of inquiry led to a description of the dis- 
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solution modes of King models within a tidal cutoff. 
We demonstrate that competition between the relax- 
ation, dynamical and stellar evolution time scales leads 
to a split between dissolutions dominated by relaxation 
processes and those dominated by dynamical processes. 
By sampling the relevant time scales, the boundary is 
mapped. 

We also generate a comparison of different stellar evo- 
lution codes linked to the same dynamics code and run 
against the same initial conditions, demonstrating that 
the specifics of the choice of stellar evolution recipes are 
amplified by the stellar dynamics, and impact the results 
of the simulation. 

The structure of this paper is as follows: in <j2] we 
describe AMUSE and its specific use for the dissolv- 
ing star cluster problem. This is followed by <|3] where 
the physical model, including details of the CW stellar 
evolution approximation, are deta iled. jS] contains the 
validation of AMUSE runs (jp~T| and j jC2| ), a study of 
the conseque nces of the variance in initial conditions on 
simulations ({ 4.3 ), an exploration of th e typ es of dissolu- 
tion which can disrupt a King model (f 4.4 1, and a direct 
comparison of stellar evolution codes Q4To|. Finally, [JH] 
summarizes the results and proposes future work. 

2. COMPUTATIONAL FRAMEWORK 

Historically, astrophysical simulation codes have been 
constructed by a single author or by a small group work- 
ing closely together. The typical course of the develop- 
ment begins with a simple solver for a specific physical 
problem (for example, an N-body integrator for a col- 
lisionless system) and then gradually extends to cover 
more varied physics (to continue the example, collisional 
physics or stellar evolution effects might be added). The 
code is monolithic: all physical processes are incorpo- 
rated into a single, unified program. Once written, 
changes to the code become progressively more and more 
difficult, as each physical model is tightly coupled to all 
the others. It is possible to extend the size of the de- 
velopment team beyond a few developers, but this adds 
a great deal of complexity to the management of the 
software creation process. The extension of the code be- 
comes limited not by the complexity of the physics in- 
volved, but by software development challenges. 

Despite these difficulties, a number of very successful 
monolithic codes have been developed. Among these are 
the N body serie s of codes (for a revi ew, see Aarseth 
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1999 ), Gadge t flSpringel et~aL| |2001D 
eil et al.||2000D and STARLAB (Portcg 



ell et al. 2000) and STARLAB (Portegies Zwart et al. 
200l| . Nevertheless, it is becoming clear that the limits 
of the monolithic approach are being reached. In order 
for new physics to be added to these packages, the pro- 
grammer (or team of programmers) must be an expert 
in the new physics being added, as well as in every phys- 
ical domain already present in the tightly coupled code. 
This, combined with the difficulty of modifying any ex- 
isting physics in these packages, limits the effectiveness 
of further work in the monolithic direction. 

The AMUSE philosophy is to move away from a 
general-purpose multiphysics "solver" and toward a suite 
of standardized special-purpose "evolver" modules. Each 
evolver knows about only a single physical domain, and is 
responsible for advancing a known system state through 
time by implementing the physics specific to that do- 



main. In particular, an evolver is not expected to take 
into account any physics outside its own domain in its 
calculations. 

The AMUSE standard defines four physical domains 
of interest: gravitational dynamics, stellar evolution, hy- 
drodynamics and radiative transfer. A standard interface 
to an evolver is defined for each of these domains. For 
example, the stellar dynamics interface specifies how par- 
ticles are communicated to the evolver (added, removed, 
updated) and how to make the evolver step forward a 
given number of time units. Similarly, the stellar evolu- 
tion interface specifies how to communicate star proper- 
ties (mass, age, metallicity, etc.) to the evolver and how 
to make the evolver advance to a specified time. 

All evolvers for a given physical domain are accessible 
within the AMUSE environment through this standard 
interface. This means that evolvers within a domain are 
interchangeable. As shown in <Q it is possible for a re- 
searcher to switch between several stellar evolution mod- 
els to test the effect of changing the physical approxima- 
tions used on the behaviour of the entire system. The 
same is true of the other physical domains. 

Wherever possible, the AMUSE approach is to re- 
use existing codes instead of writing new ones. This 
means that many special-purpose stand-alone solvers can 
be turned into AMUSE modules. The framework pro- 
vides a quick and easy method for wrapping an exist- 
ing code in one of the standard interfaces and making 
it available within the AMUSE environment. At the 
discretion of the author, such a module may also become 
part of the AMUSE package distributed on the web to 
interested researchers. Alternately, it is possible to cre- 
ate a "private" AMUSE module which exists only on 
the author's computers. 

In order to make use of AMUSE , the researcher writes 
a "top-level" script (using the Python scripting language) 
which instantiates a set of evolvers relevant to the prob- 
lem being studied. All communication and synchronisa- 
tion between the evolvers are handled by this script. In 
this work, the top-level script creates a stellar dynamics 
evolver (in our case, an N-body code) and a stellar evo- 
lution evolver. It then begins a loop in which dynamics 
and evolution are advanced in tandem, with synchroniza- 
tion between them as needed. It also implements a tidal 
cut-off by removing escapers from the simulation at fixed 
time intervals. 

AMUSE uses the messa ge-passing interface MPI (see, 
for example, Walker 1994]) to allow each evolver to exist 
in its own process, possibly in parallel and on a different 
machine than the controlling Python script. Each evolver 
is written in the language of choice of its original author. 
Already present in AMUSE are modules written in C, 
C++, Python, Fortran and Java. MPI was chosen based 
on the experience of MUSE (which used Swig and f2py 
instead), and allows for both parallelization and for each 
module to reside in its own, unique namespace. AMUSE 
is capable of running on a gr id for massively parallel cal- 
culations ( |Drost et aL]|2012[ ) 

Table IT] lists the specific AMUSE modules used in 
this work? The PH4 evolver provides N-body dynamics 
using Sapporo, a GRAPE emulator, for GPU accelera- 
tion. ph4 is an MPI-parallel, adaptive block time-step, 
GRAPE-accelerated, 4 th -order Her mite integrator, sim - 
ilar in construction to PHiGRAPE ( |Harfst et aL]|2007[ ). 
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Table 1 

AMUSE Modules Used In This Work 


Module 


Type 


Reference 


ph4 


N-Body Dynamics 


McMillan ct al. (2011) 


Sapporo 


GRAPE Emulator on GPU 


Gaburov ct al. (2009) 


SSE 


Stellar Evolution 


Hurley ct al. (2000) 


EFT89 


Stellar Evolution 


Egglcton ct al. (1989) 


VSSE 


Stellar Evolution 


Chcrnoff & Weinberg (1990) 


SeBa 


Stellar Evolution 


Portcgics Zwart & Vcrbunt ( 1996 ) 



SS E provides stellar evolution. knnCUDA (Garcia et 



al. 2008 Garcia 20081 is used to compute densities (us 
ing 12 tn nearest neig hbours) in a standalone c ode similar 
to that described by Casertano & Hut (19851, but sepa- 
rate from AMUSE . This code uses an exact algorithm 
to find all nearest neighbours, regardless of distance. 

Additionally, the EFT89 module, the SeBa module 
and the VSSE (Very Simple Stellar Evolution) module 
were used to provide the simplified stellar evolution mod- 
els explored in section [4] VSSE was written specifically 
for this work, and is designed to allow a researcher to eas- 
ily add simple analytical stellar evoluti on m odels, here 
the prescription of CW (see also section 3.3). 



The cluster was allowed to evolve dynamically using 
the adaptive, block time-step algorithm embedded in 
ph4. At regular intervals of 1 Myr (chosen to resolve 
mass-loss processes in t he most massive stars, see |Porte^1 
gies Zwart et al. (1999)), the stellar evolution mass loss 
was computed and applied to the synchronized N-body 
model. Various initial conditions have differing dynami- 
cal times, so the number of synchronizations per dynam- 
ical time varies between them. 

The specific software architecture used in this work is 
sketched in Figure[l] Note that GPU acceleration is used 
for both N-body dynamics and for nearest neighbour cal- 
culation. 

Our runs were conducted on a cluster of 24 dual-socket 
Intel Xeon X5650, 2.66 GHz nodes with a total of 12 
cores each. Each node contains 64 GB of RAM and six 
nVidia GTX 480 GPUs. Our runs ran on a single node, 
using two GPUs: one for dynamics calculation and one 
for density estimation. A typical run would use up to 
three cores: one each for process control, dynamics su- 
pervision, and stellar evolution. 

3. PHYSICAL MODEL 

3.1. Initial Conditions 



King (1966) models with W = 3 and W — 7 were 
used as initial conditions for our model clusters, repre- 
senting relatively diffuse and relatively centrally concen- 
trated systems. A simple power-law stellar mass function 



dN oc m a dm 



(1) 



was used in conjunction with a random number generator 
to assign masses to each particle. Following CW, the 
slope of the mass function was taken to be either a = 1.5 
or a — 2.5 with masses ranging from O.4M to 15M Q . 
A Salpeter mass function would correspond to a = 2.35. 

Runs are grouped by "family," a parameter also de- 
fined by Chernoff & Weinberg. The family parameter 
fixes the relaxation time at the tidal radius, which ef- 
fectively changes the ratio of the stellar evolution and 



dynamical time scales for a given model. Four values of 
this relaxation time are chosen, as summarized in Table 
[2j The relaxation time within the tidal radius is 



t 



rlx,CW 



G^m* In AT 



(2.57 Myr) F t 



cw 



(2) 



where r t is the King truncation radius and m* is a typ- 
ical stellar mass (here taken to be M Q ). By making the 
assumption that r t is equal to the Jacobi radius rj of the 
cluster, Few is defined to be 



Ft 



cw 



M R g 220 km s" 1 1 



M Q kpc 



InV 



(3) 



where R g and v g are the radius and velocity of the clus- 
ter, assuming a circular orbit about a parent galaxy (see 
CW). In generating our initial conditions, the choice 
v g = 220 km s _1 is made and R g is computed to match 
one of the specified families. 
The more conventional half-relaxation time is given by 

M l/2 r 3/2 

^-0.138 G1/2 (m)lniV 

= (3.54xl0 5 yr)F cw f^V /2 ff (4) 
\ r t J (m) 

where (m) is the average mass of a star in the cluster 
(TPZ). For a = 1.5 or 2.5, (m) = 2.54M or 1.01M o 
respectively. The properties of each family are listed in 
Table [2j For convenience, a brief summary of the relevant 
King model parameters is also provided in Table [3j 

The initial conditions were generated using Starlab 
, tools and are stored in the Starlab format, which is 
AMUSE -compatible. 

3.2. Tidal Truncation 

Tides are simulated by simple truncation at the Jacobi 
radius 



r.j 



M 



3AL 



alaxy 



1/3 



Ro 



(5) 



which is assumed to be equal to the King truncation ra- 
dius rt- Once every dynamical time any stars which pass 
beyond the Jacobi radius are removed from the simula- 
tion. The dynamical time for the cluster is 



t-dyn 



GM b / 2 



(-2U) 



3/2 



(6) 



wher e U is the gravitational pote ntial energy of the sys- 
tem (Portcgie s Zwart et al.||2010| . 
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Density Interface 
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GPU 



Figure 1. The AMUSE software architecture, showing the particular modules employed for the bulk of this work and the hardware 
acceleration used. 



Table 2 

Family Summary 











trh (Gyr) 






Family 


Few 


trlx.CW 




a = 1.5 




a = 2.5 






(Gyr) 


W =3 


W = 7 


W =3 


Wo = 7 


1 


5.00 x 10 4 


128 


1.01 


0.300 


2.46 


0.728 


2 


1.32 x 10 5 


339 


2.68 


0.793 


6.49 


1.92 


3 


2.25 x 10 5 


577 


4.56 


1.35 


11.1 


3.28 


4 


5.93 x 10 5 


1522 


12.0 


3.56 


29.2 


8.64 



Table 3 

Relevant King Model Parameters 



W 


c 




p(oy (p) 


3 


0.67 


0.27 


84 


7 


1.53 


0.12 


6430 



3.3. Very Simple Stellar Evolution (VSSE) 

For direct comparison with CW and TPZ, we imple- 
mented the ChernofF & Weinberg stellar evolution model 
in AMUSE . In this model, a main sequence star re- 
mains at constant mass until it has exceeded its lifetime. 
At that point, the star is transformed into a remnant of 
lower mass, and evolution ceases. 

The stellar lifetime is determined by fitting a cubic 
spline to the data points listed in Tabl e |4l which were 
taken by Chernoff & Weinberg from Miller fc Scalo 



(1979). The remna nt mass is derived usi ng Equation 7[ 



which was based on Iben & Renzini ( 1983 ) . Note that the 
lifetimes listed are assumed to correspond to Population 
I stars. 



'O.58M +0.22 {rm-Mo). 
0.0, 

1.4M , 



Table 4 

VSSE Lifetimes 



Initial Mass 
(log 10 [m;/M Q ]) 


Main Sequence Lifetime 
(logio [tse/yr]) 


1.79 


6.50 


1.55 


6.57 


1.33 


6.76 


1.11 


7.02 


0.91 


7.33 


0.72 


7.68 


0.54 


8.11 


0.40 


8.50 


0.27 


8.90 


0.16 


9.28 


0.07 


9.63 


-0.01 


9.93 


-0.08 


10.18 



m, < 4.7M Q 
4.7M© <m t < 8.0M o 
8.OM < m t < 15.0M Q 
(7) 



RESULTS AND DISCUSSION 
4.1. Validation 



We first set out to test whether or not the AMUSE 
framework could reproduce known results. CW produced 
a well studied catalog of Fokker-Planck code simulations. 
TPZ performed runs for the same initial parameters us- 
ing both a Fokker-Planck code and an N-body code, cal- 
ibrating the removal of escapers in the former against 
the latter. Our initial goal was to see if, and how well, 
AMUSE could replicate the TPZ N-body results. 

In TPZ, curves for four choices of parameters are pub- 
lished. Using the notation (Wo, a, family), they are: (3, 
2.5, 1), (3, 2.5, 4), (7, 1.5, 4) and (7, 2.5, 1). We there- 
fore generated initial conditions for each of these cases 
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Figure 2. Evolution of the mass of a cluster of N=32,000 stars, with initial cluster parameters (Wq, a, family), from left to right: (3, 
2.5, 1), (3, 2.5, 4), (7,1.5, 4) and ( 7, 2.5, 1). The track produced by A MUSE is shown as the solid line (blue in the electronic version). 
The corresponding run result from |Takahashi fc Portegies Z wart (200u| (N-Body model) is shown as the dashed line (red in the electronic 
version). 



and ran them using AMUSE . 

Figure [2 shows the result of our attempt to reproduce 
selected TPZ runs. The qualitative shape of the curves 
are similar, but it is difficult to make any quantitative 
statement about the agreement or disagreement between 
the two sets of simulations. 

An obvious suspect for these variances is the difference 
in stellar evolution recipes, which are known to be simi- 



lar but not ide ntical. TPZ used SeBa (Portegies Zwart & 



Verbunt 



f 996 ) to model stellar evolution, while we have 



used SSE. it is likely that small differences in the han- 
dling of main sequence lifetime and remnant masses for 
short-lived O and B stars could alter the later evolution 
of the cluste r due to early mass loss. This is not the only 
possibility. §4.5| reviews the variance introduced by the 
choice of stellar evolution model. 

One known difference is in the handling of kicks in 
supernova formation. Our top-level AMUSE script ig- 
nores any kick prescriptions built into the stellar evolu- 
tion codes. However, some of the models (for example, 
CW and SSE) leave a zero-mass "remnant" for some ini- 
tial stellar masses. This represents the detonation of the 
remnant during the supernova, and the ejection of any re- 



maining material. Our AMUSE script treats this as the 
star disappearing. Slight differences in these prescrip- 
tions are likely to contribute to variations seen between 
stellar evolution codes. It is unlikely, however, that this 
was the dominant source of variation. 

This led us to consider the inherent spread in our re- 
sults due to random variations in the initial conditions 
(individual star masses, positions, and velocities), as op- 
posed to systematic differences between the integrators. 
We conducted multiple random realizations of each of 
our chosen initial parameter sets. Each parameter set 
(for example, Wo — 3, a — 2.5, Family = I) describes a 
continuous model, and we transform this model into a 
discrete mass spectrum and set of spatial and velocity 
data using a variable random number seed. We refer to 
a set of random realizations of a given parameter set as 
"formally equivalent." 

Figure [3 shows the results of 21 different realizations 
of three selected models. These models are the (Wo, a) 
combinations (3, 1.5), (7, 2.5) and (3, 2.5) with all four 
families simulated for each case. This subset is 12 of the 
16 possible combinations of Wo, a and family for the 
TPZ choices of these parameters. We do not show (7, 
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Figure 3. Evolution of the mass of a cluster of N=32,000 stars, with initial cluster parameters (Wo, a), from top-left clockwise: (3, 1.5), 
(7, 2.5), and (3, 2.5). For each family, the solid line indicates the median value, the dashed lines indicated the 25 th and 75 th percentiles, 
and the shaded region indicates the total parameter space explored. 
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1.5) as it is largely uninteresting: the cluster dissolves 
quickly, much like the (3, 1.5) case. Furthermore, this is 
not a useful parameter set for validation because TPZ do 
not publish their cluster mass evolution curves for these 
combinations of parameters. 

The variation due to randomness in the initial con- 
ditions is small, except for the case Wo = 3, a = 2.5, 
where the cluster dissolves, but "slowly" relative to stel- 
lar evolution. "Slowly" refers to the fact that the ini- 
tial mass loss driven by supernovae of O and B stars is 
not sufficient to disrupt the cluster. We explore this in 
further detail in |4.4| For Wo — 3, a = 1.5, the early 
stellar evolution loss due to massive stars dominates the 
cluster's evolution, and it dissolves before differences in 
initial conditions can have much effect. Conversely, for 
Wq = 7, a = 2.5, the cluster is quite long-lived and the 
effects of initial condition variability are smoothed out 
over time. 

4.2. Comparison with TPZ 

Figure [4] also shows curves derived from TPZ's Star- 
LAB runs overplotted on an AMUSE simulation made 
using 21 runs done with ph4 and the SeBa module. De- 
spite using the "same" stellar evolution recipe, the re- 
sults clearly do not agree. This is due to "drift" in the 
SeBa code since the TPZ paper was published in 2000. 
"Drift" refers to changes in the code underlying SeBa 
over time as the stellar evolution model is improved. Fig- 
ure [5] shows one such change: the remnant mass kept by 
the code has changed in the 12 years between the publi- 
cation of TPZ and the current SeBa. 

The lesson to simulators is clear: codes can change 
with time, and these changes can produce significant dif- 
ferences in the results of simulations. 

4.3. Sources of Variance 

We would like to know why some of the curves in Fig. 
[3] have such large spreads in their tracks while others 
are constrained to a narrow area. In particular, the 
Wq — 3; a = 2.5 case behaves quite differently from the 
other cases. There are two obvious places in the dis- 
cretization process where random variations might play 
a decisive role — the masses of the most massive stars, 
which explode early and can eject a significant amount of 
material from the cluster, and their locations, since mass 
ejected from the cluster center is much more destruc- 
tive t o the cluster than mass ejected from the outer re 



gions ( Vesperini et al.|20 09; Pelupessy & Portegies Zwart 



2012 1. 

To separate these two effects, we explored the conse- 
quences of holding one of these "random" parameters 
constant and varying the other. For each of the two pa- 
rameters we ran 11 simulations with random choices of 
the other. Figure [6] illustrates the results. It is clear from 
the plot that varying the masses while holding the po- 
sitions fixed has a larger effect than holding the masses 
constant and varying the positions. The mass spectrum 
effect is larger than the spatial effect by about a factor 
of about 2. These two effects are of the same order, and 
neither is negligible when comparing simulation results. 

4.4. Types of Dissolution 



We next turned to the question of what it means for a 
star cluster to dissolve "slowly." There are three compet- 
ing time scales in this simulation: the stellar mass loss 
time scale tsE (which is of order a few tens of Myr at 
the beginning of the simulation, and scales proportion- 
ally to t as the evolution of the cluster progresses), the 
dynamical time scale td yn (defined in equation J6J) and 
the relaxation time scale t r h (defined in equation ffl. For 
a given choice of mass function slope a the stellar mass 
loss time scale is fixed; it is an intrinsic property of the 
stars themselves and not of their dynamical phase space 
configuration. 

The runs shown in Fig. [3] show two modes of cluster 
dissolution: the flat slope of the long-lived (Wo = 7, 
a = 2.5) runs or the short-lived "ski jump" curve of 
the short-lived (Wq = 3, a = 2.5) runs. The (W = 3, 
a = 1.5) runs live extremely short lives. Nevertheless, it 
is clear from examining the family 4 curve that the "ski 
jump" is present in them too, and that they are more 
similar to (Wq = 3, a = 2.5) than to (Wq = 7, a = 2.5). 
It is clear that the lifecycle of massive stars plays a role in 
the disruption of the cluster in the short-lived "ski jump" 
cases. In the long-lived case of (W — 7, a — 2.5), the 
cluster survives the early stellar mass loss and evolves 
according to relaxation processes. 

This inspired additional runs, conducted for Wo = 4, 5, 
and 6 with a = 2.5. A randomly-realized selection of 10 
runs was conducted for each family, and the results are 
plotted in Fig. J7| The case (Wo = 5, a = 2.5, Family=4) 
is particularly mteresting since it dissolves via dynamical 
processes while the remainder of the (Wo = 5, a = 2.5) 
models dissolve via relaxation processes. 

Figure [8] is a schematic of the relationship between the 
three time scales for N — 32k. We have subdivided the 
permitted region into two areas. In the right-most region, 
the stellar mass loss time scale is much smaller than the 
dynamical time scale. This corresponds to the behaviour 
seen in the (Wo = 3, a = 1.5) runs: the cluster is rapidly 
disrupted by the death of massive stars early in its life. 
In this region, clusters dissolve due to disruption. 

In the left-most region, the opposite is true. In the 
left-most region as the ratio of the dynamical to stellar 
mass loss timescales decreases so do the effects of stellar 
mass loss on the cluster structural evolution. For systems 
in this region (e.g. the family of models with (Wo = 7, 
a = 2.5)), the cluster dissolution is driven mainly by 
two-body relaxation mass loss. In our terminology, the 
cluster dissolves due to relaxation. 

There is a small forbidden region near the horizontal 
axis which is not shown due to its small size. This region 
is forbidden because t r i x .cw cannot be smaller than tdyn- 
Similar diagrams exist for other values of N and a. 

The intermediate region (shown as a line in Fig. [8]) is 
an interesting case. We have seen that for (W = 3, a — 
2.5) the cluster dissolves due to disruption. However, 



in published plots such as Fukushige & Heggie (1995), 
there are cases within this region of cluster dissolution 
due to relaxation. We suspected that this was a change 
in behaviour caused by a differing ratio of the relaxation 
time to the dynamical time. In order to control this ratio, 
we adjusted the number of stars N for the case (Wo = 3, 
a = 2.5, Family=l). 

The result is plotted in Fig. [9] In this plot, it appears 
that N = 2000 (or N = 1000) could be either a relaxation 
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Figure 4. Evolution of the mass of a cluster of N=32,000 stars, with initial cluster parameters Wo = 3, a = —2.5, family = 1 (left) 
or family = 4 (right). The range of variation in the AMUSE runs (made using ph4 and SeBa) is shown, with solid lines indicating the 
median, dashed lines the 25 th and 75 th percentiles, and the s haded area the entire envelope of the explored parameter space. A comparison 
is made to runs published by |Takahashi fc Porte gics Zwart (2000) in red. 




M ZAMS [M a ] 

Figure 5. Change in the handling of remnants within SeBa between TPZ (published 2000) and the AMUSE SeBa module. Note that 
SeBa has been improved to better match the current understanding of stellar evolution. 




t (Gyr) t (Gyr) 

Figure 6. Evolution of the mass of a cluster of N=32,000 stars, with initial cluster parameters Wo = 3, a = —2.5, family = 1. The range 
of variation seen in AMUSE runs is shown in blue with solid indicating the median, dashed lines the 25 th and 75 th percentiles, and the 
shaded area the entire envelope of explored parameter space. The left-hand plot shows the result when the same mass list is applied to 
differing spatial positions, while the right-hand plot shows the result when the same spatial positions are matched to differing mass lists. 
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Figure 7. Evolution of the mass of a cluster of N=32,000 stars, with initial cluster parameters (Wo, a), from top-left clockwise: (4, 2.5), 
(6, 2.5), and (5, 2.5). For each family, the solid line indicates the median value, the dashed lines indicated the 25*' 1 and 75 th percentiles, 
and the shaded region indicates the total parameter space explored. The families are 1, 2, 3, and 4 in order from left to right. Note that 
family 4 in the Wo = 5 case shows dynamical dissolution behaviour while the other three families show relaxation dissolution behaviour. 
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Figure 8. The space of combinations of the three fundamental star cluster time scales for N = 32k. We identify two permitted regions 
where the death of the star cluster can be characterized as due to "disruption" (cluster dissolution triggered by the death of massive stars) 
or due to "relaxation" effects (due to the slow evolution of the cluster's dynamical parameters). The label for each point is (N, Wo, -a, 
Family). Dissolution results are denoted with diamonds. Relaxation results are shown as crosses. 




Time [Gyr] 



Figure 9. Star cluster mass loss tracks for clusters with variable number of stars N. Each N has 11 random realizations. The solid 
lines show the median of the runs while the dashed lines show the 25th and 75th percentiles. All initial conditions are King Models with 
W = 3, a = 2.5, Family = 1. 
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or disruption curve. However, if it were to be a relaxation 
curve, we would expect to see a longer lifetime for the 
cluster than the N = 1000 case (comparing, for example, 
to Fig. I3J. The N = 2000 lifetime is shorter than the 
N = 1000 lifetime, though, which places the N = 2000 
case in the disruption region. A similar argument holds 
for the N = 1000 case. These small N runs are, in 
fact, more analogous to the (Wo = 3, a = 1.5) runs 
from Fig. [3j In these cases, the cluster is dissolving in 
a small number of dynamical times. This set of runs 
therefore shows that all of the dissolutions for this set of 
parameters have a dynamical disruption character. 

One should note that dynamical disruption does not 
have to be a rapid process relative to a Hubble time. 
While there are certainly cases of clusters being disrupted 
within a few tens of Myr, there are also dynamical dis- 
ruption cases where the cluster survives for several Gyr. 

4.5. Stellar Evolution Comparison 



Figure HO] shows a single N = 32, 000 model (W = 3, 
a = 2.5, family = 1) evolved using different stellar evo- 
lution models. AMUSE easily allows switching stellar 
evolution models in the same code. We used the SSE 
module, as used in the previous sections, as well as the 
idealized "VSSE" ("very simple stellar evolution") mod- 
ule described pr eviously, and an imple mentation of the 
recipes given by Eggleton et al. (1989). All curves used 
were computed by AMUSE using PH4 for gravitational 
dynamics. Also included is a plot of the population syn- 
thesis (stellar evolution only) results for these initial con- 
ditions using the different stellar evolution models. 

The small differences in the mass loss rates of the stars 
in the various models are amplified by the effects of 
gravitational dynamics and tidal stripping to produce 
differences in overall cluster lifetime of approximately 
25%. The early supernovae experienced by O and B stars 
drives the cluster toward a faster dissolution because the 
relationship between the mass loss and shrinking tidal 
radius is self-reinforcing. When mass is lost, the tidal 
radius shrinks, which may leave some stars beyond the 
"edge" of the cluster, thus leading to additional mass 
loss and starting the cycle over again. These are system- 
atic differences comparable in magnitude to the random 
run-to-run variations discussed in the previous section. 

5. CONCLUSIONS AND FUTURE WORK 

5.1. Scientific Results 

In this paper, we examined whether or not AMUSE 
could be compared to published runs, in particular those 
of TPZ that reproduce prior work in CW. We also studied 
the source of variance between formally equivalent simu- 
lation runs, and compared the effect of varying the stel- 
lar evolution model on the lifetime of a star cluster. The 
dissolution mechanism (disruption versus relaxation) for 
star clusters was explored. 

Our AMU SE runs agree with ITakahashi & Portegics 
Zwart ( |2000 ) (which agrees with |Chernoff fc Weinberg 
Q1990p ), when all variances are taken into account. The 
sources of understood variances are: the specific differ- 
ences in the random realization of the mass spectrum, 
the random realization of the initial spatial distribution 
of stars, and the difference in remnant masses produced 
by the 2000 and 2012 versions of SeBa. These small dif- 
ferences are amplified by the dynamical interactions of 



stars over the cluster's lifetime. It is not surprising that 
the remnant mass prescriptions have changed over the 
past 12 years, as this has been a topic of active research 
in the community. 

The direct comparison of stellar evolution codes has 
yielded an interesting result: prescriptions which differ 
by no more than a few percent in population synthesis 
studies can drive otherwise identical star cluster simula- 
tions to evolve at paces that differ by up to 25%. The 
cause is an amplification effect between mass loss due to 
stellar evolution and mass loss due to tidal stripping. As 
the cluster loses mass, particularly due to early super- 
novae of O and B stars, the tidal radius of the cluster 
shrinks and more outlying stars are stripped. 

The use of random realization to generate initial con- 
ditions was also examined, and we found that formally 
equivalent (i.e.: same Wq, a and CW family) initial 
conditions do not always follow the same evolutionary 
tracks. This effect is important in those cases where the 
star cluster neither dissolves rapidly (such as the loosely 
bound Wo — 3 case with a top-heavy mass function of 
slope a = 1.5) nor remains tightly bound over a Hubble 
time (such as Wq = 7; a = 2.5). Only in the intermediate 
regime (i.e.: Wq = 3; a = 2.5) does this effect matter. 

Variation in the evolution and dissolution time of a star 
cluster in the intermediate range is due primarily to the 
variation in masses of the most massive cluster members, 
but the effect of the differences in stellar spatial distri- 
bution due to random realization is also non-negligible. 
As had been expected, the population of O and B stars 
dominates the early cluster evolution because the effects 
of mass loss from supernovae are amplified as the clus- 
ter ages. The specific random realization of the spatial 
distribution also affects the overall cluster evolutionary 
track, but only at about half of the significance of the 
top-end of the mass spectrum. This leads us to believe 
that the fraction of O and B stars alone is not sufficient 
to describe the variance they introduce; we must also 
consider the position within the cluster of those stars. 

We have shown that there is a parameter-space bound- 
ary between King models that dissolve via dynamical 
processes and those that dissolve via relaxation pro- 
cesses. Figure [8] shows the location of this boundary 
in the space of the three relevant time scales: the stellar 
evolution mass loss time scale, the cluster's relaxation 
time scale and the cluster's dynamical time scale. 

5.2. Computational Methods Results 

The AMUSE framework shows promise as a new 
method of binding together domain-specific physics codes 
to form a larger simulation. Future work will be directed 
toward improving the abilities of this young framework 
relative to existing monolithic codes. Specifically, new 
dynamical modules have recently been added to handle 
close encounters between stars as well as the formation, 
evolution and destruction of binary and multiple stars. 
The addition of this module will allow AMUSE to fol- 
low the evolution of a star cluster into core collapse and 
beyond. 

In parallel to this development, the AMUSE team has 
also implemented modules to provide gas dynamics (gen- 
eralized hydrodynamics, in fact) and radiative transfer 
processes. Using these modules, AMUSE should be able 
to perform production quality simulations including all 
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Figure 10. The evolution of the mass of a cluster of N=32,000 stars, with initial cluster parameters Wo = 3, a = —2.5, family = 1 using 
different stellar evolution models. The left-hand frame shows the results if gravitational dynamics, stellar evolution and a tidal cut-off arc 
used. The right-hand plot shows the result of stellar evolution actin g in isolation. The blue, gree n, and black l ines correspond, respec tively, 
to the stellar evolution prescriptions of |Hurley et al.| l |2000[ | (SSE), jChernoff fc We inberg 1990:) (VSSE), andjEg gleton et al.| |l989l l. 



of these components. 

The modular structure of AMUSE facilitates com- 
parison of physics modules and enables exploration of 
assumptions and approximations that is difficult or im- 
possible with other simulation codes. We have used 
AMUSE to compare the effect of changing the stellar 
evolution prescription on an otherwise identical simula- 
tion, and used that to demonstrate that "small" differ- 
ences between prescriptions are, in some cases, significant 
in the cluster's evolution. The ability to directly compare 
individual scientific codes within a multi-physics simula- 
tion is novel property of a modular framework. 

The interchangeability of modules benefits: 

• the users of simulation codes, who may now select 
codes from a "menu" of choices, 

• those who interpret simulation results, who can 
now begin to understand the consequences of the 
choice of prescriptions made, and 

• authors of new codes, who may focus their code 
on a single area of physics, knowing that it may 
be linked to a multi-physics environment easily, 
and that cross-validation of their work with exist- 
ing codes can be accomplished by performing con- 
trolled experiments against established modules. 
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